W metodzie implicite zwanej "backward Euler", służącej do rozwiązywania problemów niestacjonarnych, operator różniczkowy opisujący "fizykę" modelowanego zjawiska obliczany jest w "aktualnej" chwili czasowej.
\( \frac{u_{t+1}-u_t}{dt} - \mathcal{L}({u_{t+1}}) = {f_{t+1}} \)
Tak więc w każdym kroku czasowym konieczne jest rozwiązanie równania na stan następny \( u_{t+1} \).
Prześledźmy działanie metody backward Euler w przypadku stosowania izogeometrycznej metody elementów skończonych.
Żeby opracować solwer wykorzystujący metodę implicite w izogeometrycznej metodzie elementów skończonych, musimy przetransformować sformułowanie silne w sformułowanie słabe.
Przemnażamy więc nasze równanie przez funkcje testujące
\( (u_{t+1},w) - (dt * \mathcal{L}(u_{t+1}),w) = (u_t+dt* f_{t+1},w) \)
Do aproksymacji stanu naszego systemu w danej chwili czasowej użyjemy kombinacji liniowej funkcji B-spline. W tym celu wybieramy bazę dwuwymiarowych funkcji B-spline, określając wektory węzłów w kierunku osi układu współrzędnych, na przykład dwuwymiarową bazę funkcji B-spline drugiego stopnia
\( \{B^x_{i,2}(x)B^y_{j,2}(y)\}_{i=1,...,N_x;j=1,...,N_y} u_{i,j } \)
Będą one stosowane do aproksymacji symulowanego pola skalarnego aktualnej chwili czasowej
\( u(x,y;t)\approx \sum_{i=1,...,N_x;j=1,...,N_y} u^t_{i,j } B^x_i(x)B^y_j(y) \)
Podobnie do testowania użyjemy funkcji bazowych B-spline:
\( w(x,y) \in \{ B^x_k(x)B^y_l(y) \}_{k=1,...,N_x;l=1,...,N_y} w^{k,l } \)
Zakładając że operator różniczkowy \( {\cal L } \) opisujący "fizykę" jest liniowy, nasze równanie wygląda teraz następująco:
\( \sum_{i=1,...,N_x;j=1,...,N_y} u^{t+1}_{i,j } (B^x_i(x)B^y_j(y)-dt * \mathcal{L}(B^x_i(x)B^y_j(y)),B^x_k(x)B^y_l(y))=RHS \\ \forall k=1,...,N_x; l=1,...,N_y \). Nie ustalamy tutaj konkretnego problemu, wyprowadzenia zaprezentowane tutaj dotyczą dowolnego problemu fizycznego który da się zasymulować opisaną metodą. Na przykład dla problemu transportu ciepła mamy
\( \sum_{i=1,...,N_x;j=1,...,N_y} u^{t+1}_{i,j } \\ (B^x_i(x)B^y_j(y)-dt *\left(\frac{\partial^2 (B^x_i(x)B^y_j(y))}{\partial x^2}+\frac{\partial^2 (B^x_i(x)B^y_j(y))}{\partial y^2 }\right),B^x_k(x)B^y_l(y))=RHS \\ \forall k=1,...,N_x; l=1,...,N_y \).
Dzięki sformułowaniu słabemu możemy scałkować przez części
\( \sum_{i=1,...,N_x;j=1,...,N_y} u^{t+1}_{i,j } \\ \left( (B^x_i(x)B^y_j(y),B^x_k(x)B^y_l(y)) - dt *(\frac{\partial (B^x_i(x)B^y_j(y))}{\partial x},\frac{\partial (B^x_k(x)B^y_l(y))}{\partial x}) -dt* (\frac{\partial (B^x_i(x)B^y_j(y))}{\partial y},\frac{\partial (B^x_k(x)B^y_l(y))}{\partial y }) \right)=RHS \\ \forall k=1,...,N_x; l=1,...,N_y \).
Dzięki strukturze produktu tensorowego funkcji B-spline oraz dzięki temu że pochodna w kierunku y z funkcji B-spline w kierunku osi x daje 0 (ponieważ funkcje te są stałe w kierunku osi y) i podobnie dla pochodnej w kierunku y z funkcji B-spline w kierunku osi x mamy
\( \sum_{i=1,...,N_x;j=1,...,N_y} u^{t+1}_{i,j } \\ \left( (B^x_i(x)B^y_j(y),B^x_k(x)B^y_l(y)) - dt *(\frac{\partial B^x_i(x)}{\partial x}B^y_j(y),\frac{\partial B^x_k(x)}{\partial x}B^y_l(y)) -dt* (B^x_i(x)\frac{\partial B^y_j(y)}{\partial y},B^x_k(x)\frac{\partial B^y_l(y)}{\partial y }) \right)=RHS \\ \forall k=1,...,N_x; l=1,...,N_y \).
Zauważmy że nasz układ równań do rozwiązania to
\( \left(M_x \otimes M_y-dt* S_x \otimes M_y -dt* M_x \otimes S_y\right) u^{t+1} = F(t) \)
gdzie
\( \{ M_x \otimes M_y \}_{i,j,k,l}= \int B^x_i(x)B^x_k(x)B^y_j(y)B^y_l(y) = \int B^x_i(x)B^y_j(y)B^x_k(x)B^y_l(y) = {\bf M}_{i,j,k,l } \) to macierz masowa będąca iloczynem Kroneckera dwóch jednowymiarowych macierzy masowych,
\( \{ S_x \otimes M_y\}_{i,j,k,l} = \int \frac{\partial B^x_i(x)}{\partial x}\frac{\partial B^x_k(x)}{\partial x}B^y_j(y)B^y_l(y) = \int \frac{\partial B^x_i(x)}{\partial x}B^y_j(y)\frac{\partial B^x_k(x)}{\partial x }B^y_l(y) \) to iloczyn Kronecker jednowymiarowej macierzy sztywności i jednowymiarowej macierzy masowej, oraz
\( \{ M_x \otimes S_y \}_{i,j,k,l} = \int B^x_i(x) B^x_k(x) \frac{\partial B^y_j(y)}{\partial y} \frac{\partial B^y_l(y)}{\partial y } = \int B^x_i(x)\frac{\partial B^y_j(y)}{\partial y}B^x_k(x)\frac{\partial B^y_l(y)}{\partial y } \) to iloczyn Kroneckera jednowymiarowej macierzy masowej i jednowymiarowej macierzy sztywności.
Każdą z tych macierzy potrafimy sfaktoryzować w czasie liniowym stosując algorytm solwera zmienno-kierunkowego. Jednakże faktoryzacja ich sumy w czasie liniowym nie jest już możliwa. Jest to możliwe dopiero gdy wprowadzimy schemat kroków czasowych który umożliwia rozdzielenie macierzy na podmacierze w podkrokach czasowych tak aby koszt faktoryzacji pozostawał liniowy.
Schemat Peaceman'a-Rachford'a pozwala na rozdzielenie kroku czasowego na dwa podkroki
\( \left(M_x \otimes M_y-dt* S_x \otimes M_y\right) u^{t+1/2} = F(t+1/2)+\left(dt* M_x \otimes S_y\right) u^{t} \)
\( \left(M_x \otimes M_y-dt* M_x \otimes S_y\right) u^{t+1} = F(t+1/2)+\left(dt* S_x \otimes M_y\right) u^{t+1/2 } \)
gdzie możemy scalić macierze lewej strony w jedną macierz o strukturze produktu Kroneckera, która może zostać sfaktoryzowana w czasie liniowym za pomocą algorytmu solwera zmienno-kierunkowego
\( \left(M_x-dt* S_x \right)\otimes M_y u^{t+1/2} = F(t+1/2)+\left(dt* M_x \otimes S_y\right) u^{t} \)
\( M_x \otimes \left( M_y-dt* S_y\right) u^{t+1} = F(t+1/2)+\left(dt* S_x \otimes M_y\right) u^{t+1/2 } \)
Funkcja prawej strony reprezentuje tutaj zmiany wywołane przez siłę działającą na system w czasie kroku czasowego. Jest ona sumą dwóch elementów:

1. Stanu naszego układu w poprzedniej chwili czasowej \( (u_t,w)= \sum_{i=1,...,N_x;j=1,...,N_y} u^t_{i,j } (B^x_i(x)B^y_j(y),B^x_k(x)B^y_l(y)) \) (również przemnożonego przez funkcje testującą i scałkowanego po obszarze). Oczywiście do przedstawiania stanu naszego układu w poprzednim kroku czasowym również wykorzystujemy kombinacje liniową funkcji bazowych B-spline \( u(x,y;t)=u_t=\sum_{i=1,...,N_x;j=1,...,N_y} u^t_{i,j } B^x_i(x)B^y_j(y) \)
2. Zmiany wywołane przez siłe działającą na system w czasie kroku czasowego \( (f_{t+1},w)= (f_{t+1}(x,y),B^x_k(x)B^y_l(y)) \)